Celem niniejszej analizy jest zbudowanie modelu regresji liniowej, który będzie przewidywał poziom zadowolenia z życia (Life Ladder) na podstawie zmiennych społeczno-ekonomicznych zawartych w raporcie World Happiness Report.
World Happiness Report to coroczna publikacja, która prezentuje rankingi szczęścia narodowego, bazujące na ocenach respondentów dotyczących ich własnego życia. Raport ten koreluje te subiektywne oceny z różnymi czynnikami jakości życia, takimi jak dochód, zdrowie, wsparcie społeczne, wolność wyboru, hojność oraz postrzeganie korupcji.
Dane wykorzystywane w World Happiness Report pochodzą głównie z Gallup World Poll, czyli corocznego badania opinii publicznej przeprowadzanego w ponad 160 krajach i terytoriach, obejmujących ponad 98% dorosłej populacji na świecie. Respondenci oceniają swoje życie na skali od 0 do 10, znanej jako „drabina Cantrila”.
W niniejszym projekcie koncentrujemy się na danych zawartych w raporcie wydanym w 2024 roku, a zatem najnowsze dostępne dane pochodzą z 2023 roku.
Każdy rekord w zbiorze danych jest średnią różnych wskaźników dla danego kraju w określonym roku.
data <- read.csv("World Happiness Report 2024.csv", sep = ";")
data_clean <- na.omit(data)
dim(data)
## [1] 2363 9
dim(data_clean)
## [1] 2103 9
dim(data)[1] - dim(data_clean)[1]
## [1] 260
W oryginalnych danych są 2363 rekordy. Po usunięciu wierszy z brakującymi wartościami pozostaje 2103 rekordy. Usunięto 260 wierszy.
Dane zawierają 9 zmiennych:
Country name: Kraj, dla którego dane zostały zgromadzone.
year: Rok, w którym dane zostały zebrane.
Log GDP per capita: Logarytm PKB na osobę, odzwierciedlający dobrobyt ekonomiczny i jego wpływ na szczęście.
Social support: Miernik wskazujący poziom postrzeganego wsparcia społecznego lub dostępnych sieci wsparcia dla jednostek.
Healthy life expectancy at birth: Liczba lat, które osoba może przeżyć w dobrym zdrowiu od momentu urodzenia.
Freedom to make life choices: Miernik określający, jak wolni czują się ludzie w podejmowaniu decyzji życiowych.
Generosity: Miernik odzwierciedlający poziom hojności lub darowizn charytatywnych w kraju.
Perceptions of corruption: Miernik oceniający, jak korupcyjny uważany jest rząd, co wpływa na zaufanie i satysfakcję.
Life Ladder: Główna zmienna w analizie, mierząca subiektywne zadowolenie z życia respondentów. Jest to zmienna numeryczna, której wartości są ocenami indywidualnymi respondentów na skali od 0 do 10, gdzie 10 oznacza najlepsze możliwe życie, a 0 najgorsze. Zmienna ta będzie zależna od innych zmiennych w modelu.
class_table <- data.frame(Column = names(data_clean), Class = sapply(data_clean, class))
print(class_table)
## Column Class
## Country.name Country.name character
## year year integer
## Life.Ladder Life.Ladder numeric
## Log.GDP.per.capita Log.GDP.per.capita numeric
## Social.support Social.support numeric
## Healthy.life.expectancy.at.birth Healthy.life.expectancy.at.birth numeric
## Freedom.to.make.life.choices Freedom.to.make.life.choices numeric
## Generosity Generosity numeric
## Perceptions.of.corruption Perceptions.of.corruption numeric
Wszystkie kolumny oprócz nazwy państwa są typu liczbowego.
Podczas analizy danych ustalono, że rok nie ma istotnego wpływu na poziom wskaźnika Life Ladder. W związku z tym, w dalszych analizach zmienne nie są rozdzielane względem roku.
Aby zapewnić wiarygodność analiz i porównywalność między krajami, dane zostały przefiltrowane w taki sposób, by każdy analizowany kraj miał taką samą liczbę obserwacji. Usunięto wszystkie rekordy dotyczące krajów, które nie posiadają pełnych danych dla każdego roku z zakresu 2013–2023.
# Zakres lat
lata_docelowe <- 2013:2023
# Filtrowanie danych do wybranych lat
dane_filtered <- data_clean %>%
filter(year %in% lata_docelowe)
# Wyszukiwanie krajów z danymi we wszystkich 11 latach
kraje_pelne_dane <- dane_filtered %>%
group_by(Country.name) %>%
summarise(liczba_lat = n_distinct(year)) %>%
filter(liczba_lat == length(lata_docelowe))
# Lista krajów z pełnymi danymi
kraje_z_danymi <- kraje_pelne_dane$Country.name
# Ostateczne dane tylko dla tych krajów
dane_finalne <- dane_filtered %>%
filter(Country.name %in% kraje_z_danymi)
# Wyświetlenie liczby krajów i przykładowych danych
cat("Liczba krajów z pełnymi danymi (2013–2023):", length(kraje_z_danymi), "\n")
## Liczba krajów z pełnymi danymi (2013–2023): 82
cat("Liczba pozostałych wierszy:", nrow(dane_finalne), "\n")
## Liczba pozostałych wierszy: 902
cat("Państwa z pełnymi danymi:\n", paste(kraje_z_danymi, collapse = ", "), "\n")
## Państwa z pełnymi danymi:
## Albania, Argentina, Australia, Austria, Bangladesh, Belgium, Benin, Bolivia, Bosnia and Herzegovina, Brazil, Bulgaria, Cameroon, Canada, Chile, Colombia, Congo (Brazzaville), Costa Rica, Croatia, Denmark, Dominican Republic, Ecuador, El Salvador, Estonia, Finland, France, Gabon, Georgia, Germany, Ghana, Greece, Guinea, Hungary, India, Indonesia, Iran, Ireland, Israel, Italy, Ivory Coast, Japan, Kazakhstan, Kenya, Kyrgyzstan, Latvia, Lebanon, Lithuania, Mali, Mexico, Moldova, Mongolia, Myanmar, Nepal, Netherlands, New Zealand, Nicaragua, North Macedonia, Pakistan, Peru, Philippines, Poland, Portugal, Romania, Russia, Senegal, Serbia, Slovakia, Slovenia, South Africa, South Korea, Spain, Sweden, Tanzania, Thailand, Tunisia, Türkiye, Uganda, Ukraine, United Kingdom, United States, Uruguay, Zambia, Zimbabwe
Finalnie analizy zostaną przeprowadzone na 902 wierszach.
describe(dane_finalne %>% select(-Country.name, -year))
## vars n mean sd median trimmed mad min
## Life.Ladder 1 902 5.70 1.06 5.81 5.73 1.12 2.18
## Log.GDP.per.capita 2 902 9.65 0.98 9.75 9.71 1.13 7.57
## Social.support 3 902 0.82 0.12 0.86 0.84 0.10 0.37
## Healthy.life.expectancy.at.birth 4 902 65.26 5.70 66.51 65.81 5.54 48.80
## Freedom.to.make.life.choices 5 902 0.78 0.11 0.79 0.79 0.11 0.37
## Generosity 6 902 0.00 0.17 -0.04 -0.01 0.16 -0.34
## Perceptions.of.corruption 7 902 0.74 0.18 0.79 0.77 0.13 0.15
## max range skew kurtosis se
## Life.Ladder 7.89 5.71 -0.25 -0.44 0.04
## Log.GDP.per.capita 11.68 4.11 -0.45 -0.89 0.03
## Social.support 0.99 0.62 -1.10 0.70 0.00
## Healthy.life.expectancy.at.birth 74.60 25.80 -0.75 -0.40 0.19
## Freedom.to.make.life.choices 0.96 0.59 -0.75 0.38 0.00
## Generosity 0.70 1.04 0.99 1.43 0.01
## Perceptions.of.corruption 0.98 0.83 -1.34 1.19 0.01
vars_to_plot <- names(dane_finalne)[sapply(dane_finalne, is.numeric) & names(dane_finalne) != "year"]
# Pętla dla każdej zmiennej
for (var_name in vars_to_plot) {
variable <- dane_finalne[[var_name]]
# Ustawienie liczby słupków
num_bins <- 30
binwidth <- (max(variable, na.rm = TRUE) - min(variable, na.rm = TRUE)) / num_bins
# Histogram
histogram <- ggplot(dane_finalne, aes(x = .data[[var_name]])) +
geom_histogram(binwidth = binwidth, fill = "lightblue", color = "black") +
labs(title = paste("Histogram -", var_name), x = var_name, y = "Częstość") +
theme_minimal()
# Statystyki
stats_table <- tableGrob(data.frame(
Statystyki = c("Średnia", "Odch. std.", "Minimum", "Maksimum", "Mediana", "Skośność", "Kurtoza"),
Wartość = round(c(mean(variable, na.rm = TRUE),
sd(variable, na.rm = TRUE),
min(variable, na.rm = TRUE),
max(variable, na.rm = TRUE),
median(variable, na.rm = TRUE),
skewness(variable, na.rm = TRUE),
kurtosis(variable, na.rm = TRUE)), 3)
))
# Obok siebie: wykres i tabela
grid.arrange(histogram, stats_table, ncol = 2)
}
corrplot(cor(dane_finalne %>% select_if(is.numeric)), method = "circle", type = "upper")
set.seed(123)
train <- dane_finalne %>% slice_sample(prop = 0.9)
test <- anti_join(dane_finalne, train)
set.seed(123)
# Dodaj unikalny identyfikator wiersza
dane_finalne <- dane_finalne %>%
mutate(row_id = row_number())
# Wybierz losowo 90% danych do treningu
train <- dane_finalne %>% slice_sample(prop = 0.9)
# Dane testowe to pozostałe 10%
test <- dane_finalne %>% filter(!row_id %in% train$row_id)
# (opcjonalnie) usuń pomocniczy identyfikator
train <- select(train, -row_id)
test <- select(test, -row_id)
data_cleaned = dane_finalne %>% select(-Country.name, -year)
# metoda Hellwiga
hellwig <- function(data) { # założenie: y jest w pierwszej kolumnie
cor_matrix <- cor(data, use = "complete.obs")
R0 <- cor_matrix[1, -1] # korelacja y z predyktorami
R <- cor_matrix[-1, -1] # korelacja między predyktorami
L <- 2^length(R0) - 1 # liczba kombinacji bez pustego zbioru
comb <- expand.grid(rep(list(c(T, F)), length(R0)))
best_H <- 0
best_k <- NULL
R <- abs(R)
for (i in which(rowSums(comb) > 0)) { # pomiń pusty zestaw
k <- which(unlist(comb[i, ])) # indeksy wybranych zmiennych
H <- 0
for (j in k) {
H <- H + R0[j]^2 / sum(R[j, k])
}
if (H > best_H) {
best_H <- H
best_k <- k
}
}
return(colnames(data)[best_k + 1]) # +1 bo kolumna 1 to y
}
hellwig(data=data_cleaned)
## [1] "Log.GDP.per.capita" "Social.support"
## [3] "Healthy.life.expectancy.at.birth" "Freedom.to.make.life.choices"
## [5] "Perceptions.of.corruption"
Metoda Hellwiga pokazała, że Life Ladder w największym stopniu zależy od PKB na osobę, wsparcia społecznego, przewidywanej długości życia, poziomu korupcji w kraju oraz wolności w dokonywaniu życiowych wyborów.
model <- lm(`Life.Ladder` ~ `Log.GDP.per.capita` + `Social.support` +
`Healthy.life.expectancy.at.birth` + `Freedom.to.make.life.choices` +
Generosity + `Perceptions.of.corruption`,
data = train)
summary(model)
##
## Call:
## lm(formula = Life.Ladder ~ Log.GDP.per.capita + Social.support +
## Healthy.life.expectancy.at.birth + Freedom.to.make.life.choices +
## Generosity + Perceptions.of.corruption, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.81538 -0.29642 0.01213 0.31157 1.96475
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.682288 0.327541 -8.189 1.03e-15 ***
## Log.GDP.per.capita 0.322497 0.043407 7.430 2.79e-13 ***
## Social.support 2.549378 0.235537 10.824 < 2e-16 ***
## Healthy.life.expectancy.at.birth 0.035656 0.006361 5.605 2.86e-08 ***
## Freedom.to.make.life.choices 2.008327 0.200599 10.012 < 2e-16 ***
## Generosity -0.071480 0.117044 -0.611 0.542
## Perceptions.of.corruption -0.982369 0.130170 -7.547 1.21e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5211 on 804 degrees of freedom
## Multiple R-squared: 0.762, Adjusted R-squared: 0.7602
## F-statistic: 429.1 on 6 and 804 DF, p-value: < 2.2e-16
resettest(model)
##
## RESET test
##
## data: model
## RESET = 25.95, df1 = 2, df2 = 802, p-value = 1.201e-11
model2 <- lm(`Life.Ladder` ~ `Log.GDP.per.capita` + `Social.support` +
`Healthy.life.expectancy.at.birth` +
`Perceptions.of.corruption` +`Freedom.to.make.life.choices` ,
data = train)
summary(model2)
##
## Call:
## lm(formula = Life.Ladder ~ Log.GDP.per.capita + Social.support +
## Healthy.life.expectancy.at.birth + Perceptions.of.corruption +
## Freedom.to.make.life.choices, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.81783 -0.30322 0.00941 0.31014 1.96529
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.709053 0.324470 -8.349 2.99e-16 ***
## Log.GDP.per.capita 0.327239 0.042690 7.665 5.13e-14 ***
## Social.support 2.533030 0.233920 10.829 < 2e-16 ***
## Healthy.life.expectancy.at.birth 0.035583 0.006358 5.597 2.99e-08 ***
## Perceptions.of.corruption -0.963304 0.126322 -7.626 6.84e-14 ***
## Freedom.to.make.life.choices 1.989008 0.198012 10.045 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5209 on 805 degrees of freedom
## Multiple R-squared: 0.7619, Adjusted R-squared: 0.7604
## F-statistic: 515.2 on 5 and 805 DF, p-value: < 2.2e-16
plot(model2)
Pierwszy model sprawdza istatność statystyczną wszystkich zmiennych i ich wpływ na Life Ladder. Okazuje się, że większość z nich dobrze objaśnia poziom satysfakcji życiowej. Najmniej istotna okazała się zmienna Generosity, z tego względu zrezygnujemy z tej zmiennej w drugim modelu. Pozostałe zmienne pokrywają się ze zmiennymi wybranymi metodą Hellwiga.
# Współliniowość
vif(model2)
## Log.GDP.per.capita Social.support
## 5.221811 2.231449
## Healthy.life.expectancy.at.birth Perceptions.of.corruption
## 3.875499 1.624453
## Freedom.to.make.life.choices
## 1.525798
Badamy tutaj, czy zmienne objaśniające są ze sobą silnie skorelowane. Wspólniniowość występuje raczej niska lub umiarkowana, jednak dla poziomu PKB jest wyższa, przekracza wartość 5. Z tego względu zbadamy wartość współczynnika korelacji tej zmiennej z pozostałymi zmiennymi objaśniającymi.
# Wybierz tylko dane numeryczne, bez 'Country.name' i 'year'
num_data <- dane_finalne %>%
select(-Country.name, -year ,-row_id)
# Oblicz korelacje 'Log.GDP.per.capita' z innymi zmiennymi
correlations_gdp <- cor(num_data, use = "complete.obs")[, "Log.GDP.per.capita"]
# Usuń samą korelację z sobą
correlations_gdp <- correlations_gdp[names(correlations_gdp) != "Log.GDP.per.capita"]
# Tworzymy ramkę danych z wynikami korelacji
cor_data <- data.frame(Zmienna = names(correlations_gdp), Korelacja = round(correlations_gdp, 2))
cor_data_sorted <- cor_data %>% arrange(desc(Korelacja))
print(cor_data_sorted)
## Zmienna Korelacja
## Healthy.life.expectancy.at.birth Healthy.life.expectancy.at.birth 0.86
## Life.Ladder Life.Ladder 0.78
## Social.support Social.support 0.71
## Freedom.to.make.life.choices Freedom.to.make.life.choices 0.22
## Generosity Generosity -0.07
## Perceptions.of.corruption Perceptions.of.corruption -0.41
Rzeczywiście zmienna Log.GDP.per.capita jest bardzo silnie skorelowana ze zmiennymi Healthy.life.expectancy.at.birth oraz Social.support. Zostaje więc usunięta z modelu.
model3 <- lm(`Life.Ladder` ~ `Social.support` +
`Healthy.life.expectancy.at.birth` +
`Perceptions.of.corruption` +`Freedom.to.make.life.choices` ,
data = train)
summary(model3)
##
## Call:
## lm(formula = Life.Ladder ~ Social.support + Healthy.life.expectancy.at.birth +
## Perceptions.of.corruption + Freedom.to.make.life.choices,
## data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.72104 -0.28920 0.01071 0.29868 2.04274
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.079491 0.324956 -6.399 2.65e-10 ***
## Social.support 3.373703 0.213893 15.773 < 2e-16 ***
## Healthy.life.expectancy.at.birth 0.071612 0.004432 16.159 < 2e-16 ***
## Perceptions.of.corruption -1.272163 0.123939 -10.264 < 2e-16 ***
## Freedom.to.make.life.choices 1.618592 0.198786 8.142 1.47e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5392 on 806 degrees of freedom
## Multiple R-squared: 0.7445, Adjusted R-squared: 0.7433
## F-statistic: 587.2 on 4 and 806 DF, p-value: < 2.2e-16
plot(model3)
vif(model3)
## Social.support Healthy.life.expectancy.at.birth
## 1.740953 1.757283
## Perceptions.of.corruption Freedom.to.make.life.choices
## 1.459181 1.434928
W tym modelu wszystkie zmienne mają już wskaźnik VIF poniżej 2.
res <- model3$residuals
hist(res)
skewness(res)
## [1] -0.1254818
kurtosis(res)
## [1] 0.6823517
# Normalność reszt
shapiro.test(res)
##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.98988, p-value = 2.22e-05
qqnorm(res); qqline(res)
Sprawdzamy reszty modelu. Histogram pokazuje, że ich rozkład przypomina normalny, jest symetryczny względem zera i nie wykazuje wartości odstających, co dobrze świadczy. Test normalności Shapiro-Wilka wskazuje jednak, że powinniśmy odrzucić hipotezę zerową o normalności rozkładu. Z drugiej strony wykres kwantyl-kwantyl wygląda całkiem przyzwoicie. Większość punktów leży blisko linii prostej - oznacza to, że reszty są w przybliżeniu normalne w środku rozkładu.Odchylenia na końcach (lewa i prawa strona) — wskazują na lekkie odstępstwa w ogonach. Nie ma dramatycznych odchyleń, żadnych ekstremalnych punktów bardzo daleko od linii. Nawet gdy reszty nie mają rozkładu idealnie normalnego, to nie jest to problem, jeśli mamy dużą próbę i wykresy pozwalają na stwierdzenie podobieństwa do rozkładu normalnego. Dodatkowo naszym celem jest przecież prognoza, więc taki rozkład reszt modelu jest w porządku. Skośność bliska zeru. Tylko lekko dłuższy ogon po lewej stronie. Kurtoza niewiele ponad trzy oznacza wyszczuplony, wysoki rozkład - dużo wartości bliskich średniej, nieco cięższe ogony niż w rozkładzie normalnym.
# Heteroskedastyczność
bptest(model3)
##
## studentized Breusch-Pagan test
##
## data: model3
## BP = 122.68, df = 4, p-value < 2.2e-16
#tu chyba badamy heteroskedastyczność jeśli chodzi o PKB ale nie wiem jakie wnioski napisać
gqtest(model3, point = 0.5, data = train)
##
## Goldfeld-Quandt test
##
## data: model3
## GQ = 0.98676, df1 = 401, df2 = 400, p-value = 0.553
## alternative hypothesis: variance increases from segment 1 to 2
plot(train$Log.GDP.per.capita, train$Life.Ladder,
xlab = "Log GDP per Capita", ylab = "Life Ladder",
main = "Zależność: PKB a Szczęście", col = "steelblue", pch = 16)
abline(lm(Life.Ladder ~ Log.GDP.per.capita, data = train), col = "red", lwd = 2)
res <- residuals(lm(Life.Ladder ~ Log.GDP.per.capita, data = train))
plot(train$Log.GDP.per.capita, res,
xlab = "Log GDP per Capita", ylab = "Reszty",
main = "Reszty względem Log GDP", pch = 16, col = "grey")
abline(h = 0, col = "red")
Wynik testu Breuscha-Pagana informuje o tym, czy w modelu występuje heteroskedastyczność, czyli zmienność wariancji reszt w zależności od wartości zmiennych objaśniających — co łamie jedno z podstawowych założeń regresji liniowej. Test ten pokazał, że w naszym modelu występuje heteroskedastyczność. Przy prognozowaniu nie jest to bardzo istotnym problemem.
# Autokorelacja
dwtest(model3)
##
## Durbin-Watson test
##
## data: model3
## DW = 1.9237, p-value = 0.1384
## alternative hypothesis: true autocorrelation is greater than 0
#stabilność modelu
#Test Ramseya
resettest(model3)
##
## RESET test
##
## data: model3
## RESET = 40.041, df1 = 2, df2 = 804, p-value < 2.2e-16
sctest(model3)
##
## M-fluctuation test
##
## data: model3
## f(efp) = 0.83485, p-value = 0.965
Wynik testu Durbin-Watsona ocenia, czy w resztach modelu występuje autokorelacja — czyli zależność reszt między sobą. DW = 1,89 oznacza istotną statystycznie autokorelację reszt.
Wynik drugi pochodzi z RESET testu (Ramsey Regression Equation Specification Error Test), który sprawdza, czy w modelu regresji liniowej brakuje ważnych zmiennych lub występuje nieliniowość nieujęta w modelu. Odrzucamy H0 na korzyść hipotezy alternatywna (H₁): model jest źle wyspecyfikowany — być może brakuje jakichś istotnych zmiennych, relacje między zmiennymi są nieliniowe, ale model zakłada liniowość.
Wynik testu M-fluctuation dotyczy stabilności parametrów modelu regresji w czasie lub w innej kolejności danych (np. wg roku, kraju itp.). Model model2 jest parametrycznie stabilny — nie widać, by jego współczynniki zmieniały się istotnie w czasie lub między grupami danych. To dobry znak, szczególnie przy danych wieloletnich, bo sugeruje, że relacje między zmiennymi są trwałe i spójne.
pred <- predict(model3, newdata = test, interval="p")
mae <- mean(abs(pred - test$`Life.Ladder`))
rmse <- sqrt(mean((pred - test$`Life.Ladder`)^2))
mae
## [1] 0.8546363
rmse
## [1] 1.019808
df <- data.frame(
Actual = test$Life.Ladder,
Prediction = pred[, "fit"],
Lower = pred[, "lwr"],
Upper = pred[, "upr"]
)
df %>%
mutate(hit = if_else(Actual >= Lower & Actual <= Upper, TRUE, FALSE)) %>%
summarise(hitrate = mean(hit))
## hitrate
## 1 0.956044
MAE (Mean Absolute Error) = 0.4267 → Średni błąd bezwzględny prognozy: model myli się średnio o 0.43 jednostki w przewidywaniu poziomu szczęścia (Life.Ladder).
RMSE (Root Mean Square Error) = 0.5551 → Błąd średniokwadratowy: uwzględnia większe kary dla dużych błędów. Zarówno MAE < 0.5, jak i RMSE < 0.6, co przy skali Life.Ladder (zwykle 0–10) oznacza, że model radzi sobie dobrze z predykcją.
Różnica między RMSE a MAE jest umiarkowana, co oznacza, że nie ma dużych ekstremalnych błędów.
93% predykcji mieści się w przedziałach ufności. To naprawdę dobry wynik przy poziomie istotności 95%.
# Predykcja vs rzeczywiste
plot(test$Life.Ladder, pred[,"fit"],
xlab = "Rzeczywiste wartości szczęścia",
ylab = "Przewidywane wartości szczęścia",
main = "Porównanie rzeczywistych vs przewidywanych wartości szczęścia",
col = "blue", pch = 16)
abline(0, 1, col = "red", lwd = 2)
# Wykres reszt
ggplot(data.frame(resid = residuals(model3), fitted = fitted(model2)),
aes(x = fitted, y = resid)) +
geom_point() +
geom_hline(yintercept = 0, color = "red") +
labs(title = "Reszty vs dopasowane wartości", x = "Dopasowane", y = "Reszty")
model_year <- lm(Life.Ladder ~ Log.GDP.per.capita + Social.support +
Healthy.life.expectancy.at.birth + Perceptions.of.corruption +
factor(year), data = train)
summary(model_year)
##
## Call:
## lm(formula = Life.Ladder ~ Log.GDP.per.capita + Social.support +
## Healthy.life.expectancy.at.birth + Perceptions.of.corruption +
## factor(year), data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.22649 -0.31189 0.02231 0.32743 1.91521
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.821926 0.285369 -2.880 0.00408 **
## Log.GDP.per.capita 0.226789 0.044129 5.139 3.47e-07 ***
## Social.support 3.267151 0.238945 13.673 < 2e-16 ***
## Healthy.life.expectancy.at.birth 0.042645 0.006775 6.295 5.09e-10 ***
## Perceptions.of.corruption -1.604959 0.115291 -13.921 < 2e-16 ***
## factor(year)2014 -0.021248 0.091286 -0.233 0.81600
## factor(year)2015 0.047407 0.091668 0.517 0.60519
## factor(year)2016 -0.061083 0.091628 -0.667 0.50519
## factor(year)2017 0.029835 0.092085 0.324 0.74602
## factor(year)2018 0.106752 0.091941 1.161 0.24595
## factor(year)2019 0.096935 0.092418 1.049 0.29456
## factor(year)2020 0.089242 0.091499 0.975 0.32969
## factor(year)2021 0.055439 0.091926 0.603 0.54663
## factor(year)2022 0.011764 0.092790 0.127 0.89914
## factor(year)2023 0.108911 0.092071 1.183 0.23720
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5531 on 796 degrees of freedom
## Multiple R-squared: 0.7346, Adjusted R-squared: 0.7299
## F-statistic: 157.4 on 14 and 796 DF, p-value: < 2.2e-16
anova(model2, model_year)
## Analysis of Variance Table
##
## Model 1: Life.Ladder ~ Log.GDP.per.capita + Social.support + Healthy.life.expectancy.at.birth +
## Perceptions.of.corruption + Freedom.to.make.life.choices
## Model 2: Life.Ladder ~ Log.GDP.per.capita + Social.support + Healthy.life.expectancy.at.birth +
## Perceptions.of.corruption + factor(year)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 805 218.43
## 2 796 243.50 9 -25.067
AIC(model2, model_year)
## df AIC
## model2 7 1251.657
## model_year 16 1357.763
vif(model_year)
## GVIF Df GVIF^(1/(2*Df))
## Log.GDP.per.capita 4.949310 1 2.224705
## Social.support 2.065306 1 1.437117
## Healthy.life.expectancy.at.birth 3.903559 1 1.975743
## Perceptions.of.corruption 1.200257 1 1.095563
## factor(year) 1.052572 10 1.002565
pred <- predict(model_year, newdata = test)
mae <- mean(abs(pred - test$`Life.Ladder`))
rmse <- sqrt(mean((pred - test$`Life.Ladder`)^2))
mae
## [1] 0.4325921
rmse
## [1] 0.5405193
error <- test$Life.Ladder - pred
MAE <- mean(abs(error))
MAPE <- mean(abs(error / test$Life.Ladder))*100
RMSE <- sqrt(mean(error^2))
cat("MAE:", MAE, "\nMAPE:", MAPE, "\nRMSE:", RMSE, "\n")
## MAE: 0.4325921
## MAPE: 8.051326
## RMSE: 0.5405193
resettest(model_year)
##
## RESET test
##
## data: model_year
## RESET = 20.643, df1 = 2, df2 = 794, p-value = 1.82e-09
Próba stworzenia modelu uwzględniającego czas. Żaden rok nie jest istotny statystycznie — wszystkie mają p-value znacznie powyżej 0.05. To oznacza, że po uwzględnieniu innych zmiennych (gospodarka, zdrowie itd.), sam rok nie wnosi istotnej zmiany w poziomie szczęścia. Możliwe, że trendy czasowe są już “wychwycone” przez inne zmienne. Anova daje nam jednak informację, że RSS (Residual Sum of Squares) spadło z 545.49 do 532.57, co oznacza, że model 2 lepiej dopasowuje się do danych. Test F sprawdza, czy ta poprawa jest istotna statystycznie. p-value = 0.0004378 → oznacza, że poprawa dopasowania jest istotna na poziomie istotności 0.001. Choć żaden pojedynczy rok w modelu nie był istotny, to rok jako grupa zmiennych kategorycznych wnosi statystycznie istotną informację. Dodanie factor(year) poprawia dopasowanie modelu. Rok wnosi coś istotnego statystycznie, bo nawet niewielka poprawa dopasowania (spadek RSS) przy tak dużej liczbie danych oznacza, że efekt roku jest realny, a nie losowy.
AIC (Akaike Information Criterion) jest również niższy dla model_year. Mimo że model_year ma więcej parametrów (25 vs 7), to i tak wyraźnie lepiej dopasowuje się do danych po uwzględnieniu złożoności. To kolejny dowód (po ANOVA), że dodanie factor(year) ma sens i jest uzasadnione statystycznie.
Mimo że model_year ma lepsze dopasowanie do danych treningowych (niższy AIC, lepszy ANOVA), jego trafność predykcji na danych testowych jest minimalnie gorsza. Różnice są jednak bardzo niewielkie i mogą być przypadkowe (w granicach błędu). Dodanie factor(year) poprawia dopasowanie i ujawnia drobny efekt lat, ale nie poprawia (ani nie pogarsza znacząco) trafności prognozy.
Możemy potraktować tę część projektu jako swego rodzaju eksperyment.
data_2<-dane_finalne
data_2$Region <- countrycode(sourcevar = dane_finalne$Country.name,
origin = "country.name",
destination = "region")
data_2$RegionCode <- as.numeric(factor(data_2$Region))
levels(factor(data_2$Region))
## [1] "East Asia & Pacific" "Europe & Central Asia"
## [3] "Latin America & Caribbean" "Middle East & North Africa"
## [5] "North America" "South Asia"
## [7] "Sub-Saharan Africa"
ggplot(data_2, aes(x = factor(Region), y = Life.Ladder)) +
geom_boxplot(fill = "skyblue", alpha = 0.7, outlier.color = "red") +
stat_summary(fun = mean, geom = "point", shape = 20, size = 3, color = "darkblue",
position = position_dodge(width = 0.75)) +
labs(
title = "Rozkład szczęścia (Life Ladder) w regionach",
x = "Region",
y = "Life Ladder"
) +
theme_minimal(base_size = 14) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
data_2 %>%
group_by(year, Region) %>%
summarise(mean_life = mean(Life.Ladder, na.rm = TRUE)) %>%
ggplot(aes(x = year, y = mean_life, color = factor(Region))) +
geom_line(size = 1.2) +
geom_point(size = 2) +
labs(
title = "Średni poziom szczęścia w regionach na przestrzeni lat",
x = "Rok",
y = "Średni Life Ladder",
color = "Region"
) +
theme_minimal(base_size = 14)
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
anova_region <- aov(Life.Ladder ~ factor(Region), data = data_2)
summary(anova_region)
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Region) 6 459.8 76.63 122 <2e-16 ***
## Residuals 895 562.0 0.63
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ANOVA: p-value < 2e-16 → bardzo silne podstawy do odrzucenia H₀.
To oznacza, że region ma istotny wpływ na poziom szczęścia — przynależność regionalna wyjaśnia dużą część zmienności. Suma kwadratów dla regionu (1205) jest znacznie większa niż dla reszt (1505), co znaczy, że duża część wariancji Life Ladder jest związana z regionem. Regiony bardzo istotnie różnią się między sobą pod względem poziomu szczęścia. To idealne uzasadnienie, żeby uwzględniać region w modelach predykcyjnych lub interpretacyjnych.
Model wykazał, że największy wpływ na poziom zadowolenia z życia mają
m.in. Log GDP per capita, Social support oraz
Freedom to make life choices. Diagnoza reszt wskazuje na
pewne odstępstwa od założeń klasycznego modelu liniowego, ale model
osiąga stosunkowo niskie błędy (MAE, RMSE), co świadczy o dobrej jakości
dopasowania.
W kolejnym kroku możliwa jest rozbudowa analizy o modele panelowe z efektami stałymi lub losowymi, co może poprawić ujęcie efektów krajowych i czasowych.
```{r}}